Week 3: Basic Geospatial Operations

PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2024

Jeff Jacobs

jj1088@georgetown.edu

Wednesday, September 11, 2024

How to Do Things with Geometries

From the sf Cheatsheet

HW1 \(\rightarrow\) HW2

  • Congrats on finishing HW1! You now know how to create geometries with sf and terra
  • So now, what can you do with them?
  • For example, we’d like to be able to say things like:
    • “The new lamppost cannot be placed at \((x, y)\), since there is already a building there!”
    • “There are \(N_1\) lampposts in County 1, and \(N_2\) lampposts in County 2”
    • “The average resident in Neighborhood A lives 2 km away from their nearest bus stop

First Things First: Loading and Saving

  • Note how there were no data files in HW1 😱
  • From HW2 onwards (and in your GIS life), we’ll:
    • Download from e.g. city Open Data Portals: geo data files, but also loading on-the-fly (this week)
    • Summarize/aggregate (this week and next week)
    • Visualize findings (“Mapping Libraries” unit)

Vector Formats

Shapefiles (.shp et al.)

A shape“file” is actually (at least) three separate files bundled together:

  • Mandatory .shp: Containing feature geometries
  • Mandatory .shx: Positional indices
  • Mandatory .dbf: Data attributes
  • Optional .prj: Coordinate reference system
  • Optional .xml: Metadata

Shapefiles

Let’s see what’s inside the shapefile we first saw in Week 1, containing data on DC’s Census Tracts: Census Tracts in 2020

DC Census Tracts (with the Georgetown campus tract highlighted!) from OpenData.DC.gov

Shapefile Anatomy

From Rodrigue (2016)

GeoJSON / TopoJSON (.geojson)

  • JavaScript Object Notation: General cross-platform format
  • Useful when data is too complex for e.g. .csv
  • TopoJSON = Memory-efficient GeoJSON
  • Bonus: Inline preview on GitHub!
my_data.geojson
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [30, 20], [45, 40],
            [10, 40], [30, 20]
          ]
        ]
      },
      "properties": {
        "color": "green",
        "area": 3565747
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Polygon",
        "coordinates": [
          [
            [15, 5], [40, 10],
            [10, 20], [5, 10], 
            [15, 5]
          ]
        ]
      },
      "properties": {
        "color": "red",
        "area": 3272386
      }
    }
  ]
}

GeoPackage (.gpkg)

Raster Formats

GeoTIFF (.tif)

  • Based on TIFF format developed at NASA.

NetCDF (.nc4)

  • Used in earth sciences, as format for data sources measured and distributed multiple times per day over large full-country or full-continent areas.

Coordinate Reference Systems (CRS)

  • EPSG (European Petroleum Survey Group) Registry: Most common way to specify a CRS

    • For example, 4326 is the EPSG code for the WGS84 coordinate system

  • PROJ: Rather than opaque numeric code like EPSG, uses plaintext “proj-strings” containing parameter info: datum, ellipsoid, projection, and units (e.g. meters). Example: PROJ4 code EPSG:4326 is represented as

    +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

  • WKT: Lengthy but human-readable descriptions

Geospatial Operations

Here we’ll look at Africa:

  • The least “stretched” by different projections,
  • Interesting geometric properties like countries-inside-of-countries (South Africa, Gambia), and
  • Will be helpful for comparing/contrasting different ways of defining “neighbors”/neighborhoods

Getting the Geometries

Using rnaturalearth with mapview

France

Code 
library(rnaturalearth)
library(mapview)
france <- ne_countries(country = "France", scale = 50)
(france_map <- mapview(france, label = "geounit", legend = FALSE))

Centroid of France

Code 
library(sf)
france_cent <- sf::st_centroid(france)
france_map + mapview(france_cent, label = "Centroid", legend = FALSE)

One We Already Saw: Union

Computing the union of all geometries in the sf via sf::st_union()

Code 
library(leaflet.extras2)
africa <- ne_countries(continent = "Africa", scale = 50)
africa_union <- sf::st_union(africa)
africa_map <- mapview(africa, label="geounit", legend=FALSE)
africa_union_map <- mapview(africa_union, label="st_union(africa)", legend=FALSE)
africa_map | africa_union_map

Helpful for Rasterizing: BBox

Code 
africa_bbox <- sf::st_bbox(africa)
africa_bbox_map <- mapview(africa_bbox, label="st_bbox(africa)", legend=FALSE)
africa_map | africa_bbox_map

Convex Hulls by Country

Code 
africa_countries_cvx <- sf::st_convex_hull(africa)
africa_countries_cvx_map <- mapview(africa_countries_cvx, label="geounit", legend=FALSE)
africa_map | africa_countries_cvx_map

Convex Hull of Continent

Use st_union() first:

Code 
africa_cvx <- africa |> st_union() |> st_convex_hull()
africa_cvx_map <- mapview(africa_cvx, label="geounit", legend=FALSE)
africa_map | africa_cvx_map

One We Already Saw: Centroids

Computing the centroid of all geometries in the sf via sf::st_centroid()

Code 
africa_cents <- sf::st_centroid(africa)
africa_cents_map <- mapview(africa_cents, label="geounit", legend=FALSE)
africa_map | africa_cents_map

Pairs of Geometries

Code 
sa <- ne_countries(country = "South Africa", scale = 50)
lesotho <- ne_countries(country = "Lesotho", scale=50)
sf::st_relate(lesotho, sa)
     [,1]       
[1,] "FF2F1F212"
Code 
sf::st_covers(lesotho, sa)
Sparse geometry binary predicate list of length 1, where the predicate
was `covers'
 1: (empty)
Code 
sf::st_covered_by(lesotho, sa)
Sparse geometry binary predicate list of length 1, where the predicate
was `covered_by'
 1: (empty)
Code 
mapview(sa) + mapview(lesotho)

References

Rodrigue, Jean-Paul. 2016. The Geography of Transport Systems. Taylor & Francis.

PPOL 6805 Week 3: Geospatial Operations